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ABSTRACT 

Pseudospectral approximations in unbounded domains by Laguerre polynomials lead to 
ill-conditioned algorithms. We introduce a scaling function and appropriate numerical pro- 
cedures in order to limit these unpleasant phenomena. 
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1. Introduction 


Recently, spectral methods have been successfully applied in the approximation of dif- 
ferential boundary value problems defined in unbounded domains. At the present, different 
solution techniques are available. Among these, an approach consists in using the collocation 
method based on the nodes of Gauss formulas related to unbounded intervals. This involves 
computations with orthogonal polynomials, such as Laguerre or Hermite polynomials. For 
details about theory and numerical experiments we refer for instance to [2], [5], [6], where, 
for certain classes of problems, a spectral type convergence behavior is shown. In this paper 
we are concerned with the implementation of these methods. In fact, computations with 
Laguerre or Hermite polynomials lead to ill-conditioned algorithms and undesired round-off 
errors instabilities, even when the degree is low. After recalling some basic properties of such 
polynomials, our aim is to give suggestions in order to improve the performances in practical 
applications. 


2. Preliminary Properties and Remarks 


We review some basic properties of Laguerre polynomials. Let a > —1 be a real parameter 
and let denote the n-degree Laguerre polynomial. For any neN, is solution to the 
following Sturm-Liouville problem (see [8]): 

X lx* ^C*) + ( a + 1 “ + ni n a) ( S ) = 0 > Vle ] 0 > + °°[> ( 21 ) 

with the normalizing condition: 


L ( °\ 0) = 



( 2 . 2 ) 


_ r(n + a -fi 1) 

T(a + l)n! 

The determination of L at a given point x can be performed via the recurrence formula: 
f L& a, (x) = 1 


< 4 a) (z) = 1 + a - x (2.3) 

4 a) (x) = -[(2n + a - 1 - z)!^ 0 ^®) - (n + a - n > 2. 

n 

This procedure has a cost proportional to n. Of course, the evaluation of ^ L at a point 
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x, is given by the recurrence formula obtained by differentiating (2.3), i.e: 


jLj-.w 

dx 


Lr{x) = o 


- - 1 


d t( <x\f \ ^ T ( a ) ( \ i ~\~ OL 1 X d j(<x) { \ Tl CL Id j- (a) / \ \ O 

l Tx L ' (x) = -n iL ' (l) + n — Tx L ^ x) ' " - 2 - 

(2.4) 

Moreover, after defining the weight function w^ a \x) = x a e~ x , se]0, +oo[, a > —1, we have 
the orthogonality relation: 


f L^L$vJ a) dx = 6 n 

Jo 


r(n + a + 1) 


n! 


,n,meN. 


(2.5) 



As we shall see in the next, the expressions (2.3) and (2.4) are those mainly used in prac- 
tical applications. Unfortunately, due to the particular structure of Laguerre polynomials, 
(2.3) and (2.4) are sources of numerical instabilities. A look to the plot of the first twenty 
polynomials (see Figure 2.1) gives an idea of the troubles one can expect. In Figure 2.1, 
we have a = 0 and the plots are contained in the rectangle [0,70] x [—5000000,5000000]. 
Even for small values of the degree n, it is impossible to get a reasonable picture. Relatively 
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small values in the first part of the abscissas interval have to be compared with very sharp 
oscillations in the second part. From the theoretical point of view, we can give the following 
asymptotic results (see [8]). 

Proposition 2.1 - For any a > -1 and n > 2 the values of the relative maxima of . 
e -x/ 2 x (a4-i)/ 2 |^(a)( x )| J orrn an increasing sequence when x > Xq where : 

o a 2 — 1 

x 0 = 0 if a 2 < 1, x 0 = - — — if a > 1. 

In + a + 1 

Proposition 2.2 - For any a > —l and xq > 0, we can find a constant C > 0 such that: 

max \e~ m ^x^ +1 ^L^(x)\ « Cn'^n 0 ' 2 . (2.6) 

xc]xo,+oo[ 

From the previous statements one concludes that, even when multiplied by a decaying ex- 
ponential, Laguerre polynomials are not easy to keep under control. We will analyze in the 
next section more about this subject. 

Similar properties to those presented above are satisfied by the family of Hermite poly- 
nomials. According to [8], the n th degree Hermite polynomial H n is defined to be a solution 
of the differential equation: 

tf"(x) - 2 xH'Jx) + 2 nHjx) = 0, VxeR, (2.7) 

with the conditions: f 

#n(0) = (- 1 ) n/2 ^2)! if n is even> ( 2 ‘ 8 ) 

ifnisodd ' (2 - 9) 

Again we have a recurrence formula, namely: 

H 0 (x ) = 1 

< Hi(x) = 2x (2.10) 

. H n (x) = 2xH n -i(x) — 2(n — l)ff n _ 2 (x), n > 2. 

Orthogonality is also achieved in view of the next expression: 

f°° H n (x)H m (x)e~ xl dx = 6 nm ^2 n n\, n,meN. (2.11) 

J — oo 

Considerations similar to those previously discussed, concerning the implementation of for- 
mula (2.10), also hold. Actually, the two families are closely related by the following equali- 
ties: 

H n (x) = (-l) n ^2 n (n/2)lL { ~ / 1 2 /2 \x 2 ), if n is even, (2.12) 
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H„(x) = (-l)<-»/* 2 "((n - l)/2)!xlj; / j> )/I ( a ; J ), if n is odd. 


(2.13) 


3. Scaled Laguerre Functions 

As we noticed in Section 2, the determination of L^\x) for a given x, brings to an ill- 
conditioned algorithm. According to Propositions 2.1 and 2.2, we can try a first experiment. 
Instead of computing we could evaluate L^(x) = e~ x ^ 2 L^\x). This can be done 

by just modifying the first two terms of (2.3) by setting Zo°^(x) = e ~ x / 2 and Z[°^(x) = 
(1 + a — x)e _x / 2 . In this way, when x is large we can avoid overflow errors. The plots of 
L for a = 0 and n = 1,20 are given in Figure 3.1. This time they are all included in the 
rectangle [0,70] x [—1,1]. 



Figure 3.1 - Plot of 1J£\ for n = 1 , 20. 


Nevertheless, the new recurrence formula is still ill-conditioned. In fact, when x is large, 
the determination of Zq“^(x) gives rise to underflow problems. Our goal is to find a more 
suitable scaling function. In order to be really effective, this function, denoted by S n (x), has 
to actually depend on both x and n. Then, we would like to define: 

&n\ x ) = Sn(x)L^\x), ncN, xe]0, +oo[. (3.1) 
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We can try first by setting: ^(x) = (l + , n > 1. In this way, for any fixed n, S n l 

behaves like a polynomial and, for any fixed x, we have lim^oo S n {x) = e~ x ' 2 . After such 
a scaling, computations seem to be more flexible. An appropriate recursive formula can be 
written for Unfortunately, one can check that, with this choice of S n , the cost for 

implementing the formula is proportional to n . 

Therefore we suggest another scaling function, i.e. 

s “ ( * )= [( n r)M i+ s)] '■ "- 1 ' (s - 2) 

and we take So(x) = 1. Thus we obtain the following formula: 

' L£\x) = 1 

H*) (x ) - 4 ( a + 1 ~ x ) 

< 1 [X) ~ (a + l)(x + 4) 

[(2 n -ha - 1 - x)i2i(*) - ^ 2 - 

(3.3) 

From (3.3) we can obtain the value of L^(x) at a given x, with a cost proportional to n. In 
the next, we shall refer to L^\ obtained by (3.1) with the help of (3.2), as scaled Laguerre 
Junctions. 



Figure 3.2 - Plot of Vg\ for n = 1,20. 


£'“>(*) = 


(n + a)(4n + x) 
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It is clear that these are not polynomials. Their structure is more appealing for numerical 
tests. The plot of the first twenty scaled Laguerre functions (with a = 0) is given in Figure 
3.2. Now all the plots are included in the rectangle [0,70] x [—500,500]. The computational 
range is reduced. We make clear this fact by a numerical experiment. In the first column 
of Table 3.1 we give the values of (xJ/Zt^q (x)| for various x. There is a difference of 
13 digits between the highest and the lowest value. In the second column of Table 3.2, we 
report |Z.^(x)/i^(x)|. Now the width is reduced to 5 digits. 



X 

Column 1 

Column 2 

10 

7.273E-00 

3.129E + 03 

20 

1.779E-00 

1.035E + 05 

30 

1.607E-03 

6.215E + 03 

40 

2.405E-05 

3.743E + 03 

50 

1.954E-07 

8.379E + 02 

60 

1.786E-09 

1.554E + 02 

70 

1.638E-11 

2.285E + 01 

80 

2.187E-13 

4.015E + 00 


Table 3.1 - Comparison between Laguerre polynomials and scaled Laguerre functions. 


We remark that Z/£*)(0) = 1, VneiV. For convenience we also give the formula to compute 
the derivatives, i.e.: 


±i£\x) 


0 


££<->(*) 


4 q+5 

(*+4) 2 a+1 


fr+SRta+x) [( 2n + Q! ~ 1 

■ 4(n— l) 2 f 2(4n+x— 2) 

' 4n+x— 4 V. (4n+x)(4n+x— 4) 


- *)££&(*) - s 2g=»ifc , >(*)+ 

ii-j(x) - sifc’iM)] , ™ > 2. 


(3.4) 


We will see in the next section how to apply the scaled functions to pseudospectral compu- 
tations. 
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4. Numerical Integration 

Let N > 1 be an integer. Let us define rjfy, k - 1, . . . , N to be zeroes of L^\a > -1 
It can be shown that > 0, A: = 1, . . . , N. Then define the weights: 

T(N + a + l) V Q 


(a) 

™fc,n = 


,* = 1 JV. 


(N + 1)\[(N + 1)L^ +1 ( V { ^)] 2 
The following quadrature formula is known (see for instance [8]): 


(4.1) 


y+oo 

JO 


fw^dx = £ - / 1W) (o. 


(4.2) 


k— 1 


where 0 < t < +oo. In particular, the second term on the right hand side vanishes when 

fa} 

/ is a polynomial of degree at most 2N - 1. For convenience, we shall set r) k = r)l jj and 
w k = wjuv- We are concerned with the numerical evaluation of T) k and w k ,k = 1, ...N. An 
approximation of the zeroes of L$ can be obtained by the following procedure (see [7] and 
[8]). We first define y ki k = 1, . . . , N to be solutions of the equations. 

N — k + 3/4 


l ik ~ sin j/jt = 2?r- 


2N + a + 1 ’ 


k = 1 


Afterwards we set z k = [cos(|yjfc)] 2 , k — 1, . . . , N, and finally we get: 


Tj k ~ 2(2 N + a. + H — ■ 


4(1 - z k y i - z fc 


- 1 + 3a 2 


(4.3) 


(4.4) 


6(2iV + a + 1) 

Starting from this approximated value for 77*, , we can refine it by very few iterations of the 
Newton method. The computed zeroes are in general very accurate. However, evaluating 
4?0?i?jv) and £ L N\v^N)> k = !>•••» W, in order to use the Newton method, may give 
problems when performed by (2.3) and (2.4). Nevertheless, since the zeroes of are the 
same than the zeroes of 1$, we can use (3.3) and (3.4). This allows the determination of the 
zeroes for larger values of N . Concerning the weights, they are decreasing and converging 
to zero very fast. It is convenient to define other weights as follows (recall (3.2)): 

Wk 


w ki h = w k = 


[SNivkW 

Y(N + a + 1 )Vk 


(N + l)!T 2 (a + 1) 


4(AT+1) 


(4.5) 


_{N + a + 1)(4 N + 4 + Vk)L { MUvk)_ 

The new weights are in the machine range for larger values of N. Of course, one has to 
write: 

TV 7 U 1 1 1 > AT 1 ^ / AT _L ^ 1 / _L 1 \ \\ 

(4.6) 


T(N + a + 1) 


(N+ 1)! 


N + a f N -\- a — 1 / 


N + 1\ N V 

V 2 ) )) 
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obtaining a more reliable algorithm. 

If p is a polynomial of degree at most 2 N — 1, from (4.2) and (4.5) we obtain: 


f pw^dx = Yllp( r li') S N(Vk)]wk 

Jo k=i 


Moreover, if is a polynomial of degree at most N, we have: 


(4.7) 


f u 2 N w^ a) dx « X)[( u JV‘ S, Af)( 7 /*)] 2t£r * 
J ° k = 1 


(4.8) 


This suggests to work with u^Sn instead than when approximating the solution of a 
differential equation by the collocation method based on Laguerre polynomials. We analyze 
in the next section in which way this can be actually done. 


5. Pseudospectral Laguerre Approximations 


Let us focus our attention on a very simple equation, namely: 

— U xx + A U = / in ]0 + oo[, A > 0, 

< U{ 0) = 7eR, (5.1) 

lim*.^ = 0- 

Under suitable assumptions on / and A, one can prove that there exists a unique solution of 
(5.1) (see [2]). For approximations with Laguerre polynomials we require for U an exponential 
decay at infinity. By setting V = Ue x , problem (5.1) becomes: 

r —V xx + 2V X + (A - 1)1 7 = g in ]0,+oo[, 

' (5 2) 

l ^(0) = 1 , 


where g = fe x . 

Keeping the same notations of the previous section, let rjk, k = 1 , . . . , N be the zeroes of 
L $ . Then the solution of (5.2) is approximated by a polynomial of degree less or equal 
N satisfying the collocation problem: 


{ -vn,xx (j?*) + 2v N>x {rik) + (A - l)v N (r] k ) = g(r} k ), k = l,...,N, 
v N (0) = 7. 


(5.3) 


As usual (see [1] and [4]), (5.3) is equivalent to an appropriate linear system whose unknowns 
are %(j/h), k = 1 , . . . , N. When N tends to infinity, v n converges to V in a suitable weighted 
Sobolev space (see [2]). 
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The basic tool for computations is the derivative operator Dn in the space of polynomials 
of degree N. With the aim of getting the entries = 0, . . . , N of D n, we consider the 

Lagrange polynomials with respect to the nodes r) k , k = 0, . . . , N where we defined t ] 0 = 0). 
These are given by: 


lj(x) = 


xL^\x) 1 


l 0 (x) = 


Ln\ x ) 

ijf’(O) 


The Ij’s are polynomials of degree N . Then we have 


j = 


(5.4) 

(5.5) 


dij = i,j = 0,...,N. 


(5.6) 


Moreover for any polynomial p of degree at most iV, one gets: 


N 


^2 d ijP(Vj) = p'(Vi)> i = 0,. . . ,N. 
i— o 

The following useful relations can be recovered from (2.1): 
d 2 tM, \ rij — a — 1 d ( a ), . . 

Ly{vj) = : — L N\Vj), j = i» 


dx 2 


d 


Vi 


i 1 " (°) = -5TT iS? '(°)' 


dx 




N 


r(«)/ 




dx 2 


N(N-l) M , 
(a + !)(<* + 2) 


(5.7) 


(5.8) 

(5.9) 
(5.10) 


—L {a] (v) - 
dx * L N (7?j) ” 


'(Vj - a-2)(Vj -a + 1) iV-1 




Vi 


j = N. (5.11) 
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Thus, with the help of (5.8) - (5.11), one gets: 

Vj£L^\ Vj )Vi-Vi 
1 - a + rjj 

2r]i 


dij — 


4 a) (o) 

IjiLWiVi) 

N 

a + 1 


i,j = 




i = = 0, 


j = = 0, 


1=^ = 0. 


The second derivative operator is obtained either by squaring Dpj either by evaluating I”fa). 
In the last case, recalling (5.8) - (5.11), we have: 

[ srjy’fa) (l-a + i)i)0K-tt)-2tt .. , N ... 


fa - <*) 2 jy - 1 

37?? 3 7?, 


?«(Y 'I— J a + 1 — Vi S^l^fa) 

5 ZgJTiT 


2(jV + q+l) j#>(0) 


i = ; = l,...,iV, 


i = 1,. . . ,N,j = 0, 


j = 1, . . . ,N,i = 0, 


i _,_ 0 

(a+l)(« + 2) t-J-u. 

After taking into account the boundary conditions, we end up with the matrix corresponding 
to the linear system (5.3). 

As pointed out in the previous section, Laguerre polynomials are not suitable for com- 
putations. However we can use scaled Laguerre functions. Therefore, we define v * = 
v.Arfa)5jvfa ), k = 0 ,N. The v^s will be the new unknowns. Besides, we define a new 

A A 

matrix = {d^} as follows: 


4 = i,}= 0 N. 

b N{Vj ) 


(5.14) 
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(5.15) 


Due to (5.7) we have: 

N N A 

SNtOift'iiiVi) = £ d aMVj) = £ 1 = °» • • • ’ N ‘ 

j = 0 j=0 

Similar relations hold for higher order derivatives. Thus, by multiplying the equations in 
(5.3) by S N {vk),k = 0,...,N we can easily obtain a new linear system in the unknowns 
{, ., j = 0, . . . , N, involving the knowledge of the coefficients d %j instead than d x] . The d XJ s are 
less affected by ill-conditioning problems, since they are related to scaled Laguerre functions. 
In fact, we observe that: 

(5.16) 

(5.17) 

(5.18) 


“ l L « + 


Moreover: 


Finally, we define: 


Sn-^^n 1 ) 


N X 

(*) = E 




Q ' = )(Vi) + £ 4 k + Ci' 


i = 0,...N. 


Therefore, by substituting in (5.12), one gets. 

f OiQi 1 
VjQj V* ~ Vj 

1 - a + rii 


<kj = \ 


2vi 

Qi 

-l 

VjQj 

-N 
a + 1 


i,j = 1, • • • > Ni * ^ jt 

i =j = 1, * JV, 

i = 1,. . . , AT,j = 0, 

j = 1 = 0, 

i = j = 0. 


(5.19) 


We can argue similarly for the second derivative matrix. After solving the new preconditioned 
system we obtain i)j,j = After this, one can recover the values v n ( V] ),j 

which are in general very large, when j is close to N. On the other hand, these values are 
associated to the weight Wj that is negligible when j is close to N. Nevertheless, we remark 
that in order to evaluate the weighted norm of v N , one does not need to know its values at 
the collocation nodes. Actually, recalling (4.8) we have: 

/ + °° vjfW^dx « Yj v)™]- (5.20) 

Jo j=1 


11 


The right-hand side of (5.19) is now more suitable for computations. The procedure here 
described can be clearly generalized to other kinds of problems. 


6. Hermite Approximations 

All we developed above for Laguerre polynomials can be extended to cover the case of 
Hermite polynomials. In analogy with (2.12) and (2.13), it is convenient to define the scaled 
Hermite functions as follows: 


H n (x) — Un/l^ix 2 ), if n is even, 
H n (x) = xL^nlly^x 2 ), if n is odd. 

In particular, we obtain: 


( 6 . 1 ) 

( 6 . 2 ) 


-Hn(O) = 1, if n is even, 
H n (0) = 1, if n is odd. 


(6.3) 

(6.4) 


The derivative matrices for the Hermite case have been presented in [3]. These matrices 
can be suitably preconditioned when using scaled Hermite functions. As before, the scaling 
procedure we adopt, is only to be considered as a trick to perform computations in a better 
way. Theoretical analysis and convergence estimates remain the same, since the collocation 
scheme is not actually modified. Of course, improvements could be expected when adopting 
other functions S n in place of (3.2). 
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